#include <stdio.h>
#include <stdlib.h>
#include <math.h>


#define TIME 500	
#define OUTPATH "tauout.txt"

#define PI 3.14159

int main(){

	FILE * fp;
	fp = fopen(OUTPATH,"w");

	//float vector[] = { 0.111363,  0.015947,  0.530141,  0.163350,  0.918233,  1.228488};
	float vector[] = { 0.1,  0.03,  0.9,  0.8,  0.001,  0.01};


	float g1 	= vector[0];
	float a_1 	= vector[1];
	float a1 	= vector[2];
	float k1 	= vector[3];
	float b1	= vector[4];
	float d1	= vector[5];

	


	srand(time(0));
	
	float wsum = 0;
	int i		= 0;
	int j		= 0;

	float tau	= 0;
	int mi		= 0;
	float t	= 0;
		
	float r1;
	float r2;
	
	
	#include "model.h"
	



	float wthresh[REACTIONS];
	float wtcumsum = 0;


	//Tau fjj table
	float tmi;
	float tsigma;
	float tleap;
	float tl1, tl2;
	float teps = 1e-1;
	
	
	FILE * df;
	df = fopen("modeldif.txt", "w");
	//TODO: static
	for (i = 0; i<REACTIONS; i++)	{ for (j = 0; j < SPECIES; j++)  {dif[i][j] = -reactants[i][j]+products[i][j]; fprintf(df,"%d ", dif[i][j]);} fprintf(df, "\n");}
	fclose(df);
	
	//Constant differentials
	#include "taux.h"
	
	
	
	
	
	t = 0;
	int step = 0;
	while (t < TIME){
	
		r1 = ((float)rand())/RAND_MAX;
		r2 = ((float)rand())/RAND_MAX;
	
		//Propensities calculations
		#include "weights.cuh"
		#include "taux2.h"
	

		//Calculate wsum
		wsum=0;
		for (i = 0; i<REACTIONS; i++) {	
			//printf(" %d:%f ",i,w[i]); 
			wsum = wsum + w[i]; 
		}

		printf("wsum: %f\n",wsum);
		
		//Determine tau leap
		//TODO: popravi/static!
		tleap = 100000;
		for(i=0; i<REACTIONS; i++){
			tmi = 0; 
			for(j=0; j<REACTIONS; j++) tmi += fjj[i][j] * w[j];
			tsigma = 0; 
			for(j=0; j<REACTIONS; j++){ tsigma += fjj[i][j] * fjj[i][j]* w[j];}
			//fprintf(fp,"%f %f ",tmi, tsigma);
			//for(j=0; j<REACTIONS; j++){ fprintf(fp," %f ", fjj[i][j]); }
			//fprintf(fp,"\n");

			//printf("%f\n", tmi);
			//printf("%f\n", tsigma); 
		
			
			if(tmi > 0 || tsigma > 0){
				tl1 = teps*wsum/fabs(tmi);
				tl2 = (teps*teps*wsum*wsum)/(tsigma);
				
				if(tl2 < tl1 && tl2 < tleap) tleap = tl2;
				else if(tl1 < tl2 && tl1 < tleap) tleap = tl1;
				
				
				
				//printf("\n wsum: %f  abs(tmi): %f tsigma: %f tl1:%f tl2:%f i:%d Tleap: %f ",wsum, fabs(tmi), tsigma, tl1, tl2, i, tleap);
			}
			
		
			//tleap = (tl1 > 0 && tl1 < tleap) * tl1 + (tl1>=tleap) * tleap;
			//tleap = (tl2 > 0 && tl2 < tleap) * tl2 + (tl2>=tleap) * tleap;

		}	

		tleap = 0.007;
		
		
		
		


	




	

		if (wsum <= 0){t = TIME; printf("\n\n break \n\n"); break;}
	


		printf("%d %f %f %f\n", step, t, tleap, wsum);
		
		
		int pj;
		
		float u;
		float v;
		float mean;
		float sigma;
		float normal;
		
		fprintf(fp,"%f ", t);
		for (i=0; i<SPECIES; i++) {fprintf(fp, "(%d)%d ",i,species[i]);}
		fprintf(fp,"\n");
		
		for (j = 0; j<REACTIONS; j++)
		{
			pj = 0;
			mean = w[j]*tleap;
			sigma = mean;
			//TODO: rand() float type
			u = (float)(rand() / (1.0*RAND_MAX));
			v = (float)(rand() / (1.0*RAND_MAX));
			normal = mean + (float)(sqrt(-2*log(u))*cos(2*PI*v)*sigma);
			
			pj = round(mean + sqrt(sigma)*normal);

			
			//if(pj < 0) pj = 0;

			
			printf("\t\t  %d %f %f %f %d \n", j, mean, sigma, w[j], pj);
		
		
			for (i=0; i<SPECIES; i++){ 
				if(species[i] + pj*dif[j][i] < 0){
					species[i] = 0;
					
				}	
				else{
					species[i] = species[i] + pj*dif[j][i]; 
				}
			}
		
		} 
		
		
		
		printf("\n");
		if (step > 600)		return;
			
		
		
		
		//Print to file
		//fprintf(fp,"%f\t%d\t%f\n", t, mi,w[mi]);					
																									       //10 						     %15
		//fprintf(fp,"%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",t, species[BFP], species[mCITRINE],species[P3],species[P4],species[TALA_KRAB],species[TALB_VP16],species[iPC],w[20], mi, w[mi],species[TALA_VP16], species[TALB_VP16], r1, r2, w[3],w[13]);
	

		t = t+tleap;
		step++;

	

		//Switch mechanism;
	
		if(t >= swTime && swState > 0){
			//1 ... MCT(iEC) -> BFP(iPC)
			if(swState == 1){
				species[iPC] = INDMAX;
				species[iEC] = 0;
			}
			//2 ... BFP(iPC) -> MCT(iEC)
			else if(swState == 2){
				species[iPC] = 0;
				species[iEC] = INDMAX;
			}
			swState = 0;
		}
	




	}

	
	fclose(fp);
	return 0;
}
